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ABSTRACT 

In this paper, we present a semi-empirical model of isolated binary star for- 
mation. This model includes the effects of turbulence in the initial state of the 
gas, and has binary orbital parameters consistent with observation. Our funda- 
mental assumption is that the angular momenta of binary star systems is directly 
related to the net angular momentum induced by turbulence in parent molecular 
cloud cores. The primary results of this model are as follows, (i) A quantitative 
prediction of the initial width of the binary period distribution (o"i og p d = 1.6 — 2.1 
for a star formation efficiency in the range e* = 0.1 — 0.9). (ii) A robust, negative 
anticorrelation of binary period and mass ratio, (iii) A robust, positive correla- 
tion of binary period and eccentricity, (iv) A robust prediction that the binary 
separation of low-mass systems should be more closely separated than those of 
solar-mass or larger. These predictions are in good agreement with observations 
of PMS binary systems with periods P > 10 3 d, which account for the majority 
of all binaries. 

We conclude with a brief discussion of the implications of our results for 
observational and theoretical studies of multiple star formation. 

Subject headings: stars: formation, binaries:general, gravitation, hydrodynamics 

1. Introduction 

To fully elucidate the mechanisms underlying the origin of binary systems, it is crucial 
to compare theory with observations of the field and of nearby star-forming regions. 
However, a myriad of complex physical effects have been advanced by theorists to account 
for the formation and evolution of binary and multiple star systems, including fragmentation 
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(Hoyle 1953; Inutsuka & Miyama 1992), turbulence (Klein, Fisher, & McKee 2001; Bate, 
Bonnell, & Bromm 2002), ideal magnetohydrodynamics (MHD) (Galli et al. 2001), 
non-ideal MHD, including ambipolar diffusion (Mouschovias 1977; Boss 2000), radiative 
transfer (Boss et al. 2000), dust physics (Whitworth, Boffin, & Francis 1998), tidal 
torquing (Larson 2002), capture (Clarke & Pringle 1991), competitive accretion (Bonnell, 
Bate, Clarke, & Pringle 1997), and accretion disk processes (Laughlin & Bodenheimer 
1994). However, to date, no simulation has included all of these varied physical effects, and 
few have been run to a substantial evolutionary state. Moreover, as other authors have 
pointed out (Larson 2001), once one includes the effect of turbulence, star formation is 
inherently a stochastic process, implying that many such simulations need to be performed 
before their results can be meaningfully compared to observation. 

In this paper, we present a simplified, semi-empirical model of isolated binary star 
formation. This model includes the effects of turbulence in the initial state of the gas, and 
has binary orbital parameters consistent with observation. Our fundamental assumption 
is that the angular momenta of binary star systems is directly related to the net angular 
momentum induced by turbulence in parent molecular cloud cores. We do not require that 
the binary angular momentum be equal to that of its parent core; indeed, we shall find 
that angular momentum loss from the initial core is crucial in determining the median 
binary period, as previous authors suggested (Simon 1992; Bodenheimer 1995). While we 
do not in this work determine the mechanism underlying the loss of angular momentum, 
we find that we can describe a number of observed properties simply by specifying that the 
efficiency of transfer of mass and angular momentum from the initial gaseous core to the 
final binary are both constant. This model therefore provides an excellent framework with 
which to examine the importance of turbulence in the context of isolated star formation, 
and to explore the role which the stochastic nature of the initial turbulent state has on the 
statistical properties of binary systems. While it is certainly true that a number of other 
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physical effects other than turbulence and angular momentum dissipation may play a role 
in the star formation process, our model can elucidate their relative importance by serving 
as a strawman against which other results can be compared. 

The paper begins with an overview of observational evidence (§2), including both 
pre-stellar molecular cloud core (§2.1) and binary (§2.2) properties. Next, we describe 
our pre-stellar core model in §3. In §4, we present the methodology which connects the 
pre-stellar core model with the binary properties. The results of this model are detailed 
in §5. We conclude the main body of the paper with a discussion of the implications of 
this semi-empirical model for observational and theoretical studies of star formation (§6). 
Lastly, appendix A includes an elementary derivation of characteristic binary scales, in 
terms of their orbital parameters. 



2. Observational Evidence 
2.1. Pre-Stellar Core Properties 

2.1.1. Density Structure 

The observations of Motte & Andre (2001) consisted of a complete 1.3 mm survey of 
both pre-stellar cores and protostellar envelopes using the IRAM 30m telescope and the 
MPIfR bolometer array (MAMBO). Their observations resolved structures from 1,500 AU 
to 15,000 AU. Motte & Andre (2001) concluded that their observations of pre-stellar core 
backgrounds were flatter than predicted by a singular isothermal sphere model (Shu 1977), 
but were consistent with Bonnor-Ebert spheres (Bonnor 1957; Ebert 1955). 
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2.1.2. Velocity Structure on the Molecular Cloud Core Scale 

In 1981, in an extension of an earlier paper (Larson 1979) focusing on the large-scale 
structure of the interstellar medium, Larson (1981) published a seminal paper on the 
internal velocity dispersion in molecular clouds. Using data already published in the 
literature in a wide variety of studies, he accumulated about 50 data points for a number of 
different star forming regions, on scales ranging from 1CT 1 to 10 2 pc. He established scaling 
laws for both the mean density and the total internal velocity dispersion over projected 
lengthscales on the sky. Plotting the total internal velocity dispersion (sum of thermal and 
non-thermal linewidths) versus the maximum projected linear size of the region, he found 
that the internal velocity dispersion obeyed a power law of the form 

<7 1D (km s- 1 ) = 1.10 L(pc) - 38 (1) 

Here we have designated the one-dimensional linewidth dispersion, as inferred from 
velocities along the line of sight by o"ib, to properly distinguish it from the fully 3D velocity 
dispersion. 

Significantly, the power law nature of the turbulent linewidth-size scaling laws supports 
the premise, which Larson originally suggested, that over a wide range of scales, the 
observed interstellar turbulence is part of a scale-free hierarchy of turbulent eddies. This 
scale-free cascade should extend all the way down to about the characteristic scale where 
such eddies are damped - either the minimum Jeans mass in the cloud (Larson 1995) or 
the ion-neutral damping length (Myers & Lazarian 1998), depending on whether thermal 
or magnetic damping effects dominate, which in turn depends on the relative strength of 
the magnetic field and the ionization state of the gas. 

Subsequent authors investigated the linewidth-size scaling relationship in more detail 
for molecular cloud cores using more precise techniques. In particular, Larson's original 
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study contained only 4 cores observed with rather poor angular resolution, and in the 
absence of actual data, made blanket assumptions regarding the kinetic temperature of the 
gas used in computing the local soundspeed. Leung, Kutner, & Mead (1982) and Myers 
(1983) re-investigated Larson's linewidth- size relationship and found the exponent closer to 
1/2, rather than 1/3. In the most extensive study of molecular cloud cores to date, Jijina, 
Myers, & Adams (1999) studied 264 cores mapped in NH 3 in a wide variety of regions 
and environmental conditions. They conclude that the linewidth-size relationship for all 
cores has an exponent of .63 ± .10, though a slight difference of low statistical significance 
existed between subsamples with embedded IRAS sources (exponent .49 ± .12) and without 
(exponent .83±.18). Moreover, Jijina, Myers, & Adams (1999) found a significant variation 
in the median value of the non-thermal linewidth Av^t across regions : from a value of 0.22 
km/s in Taurus (thermal linewidth median Avt = -44 km/s), to a non-thermal linewidth 
Avnt median value of .86 km/s in Orion A (thermal linewidth median Avt = -58 km/s). 
Over all regions, cores without IRAS embedded sources or clusters had Avnt / ' Avt = -651/25 
(within one quartile of the median). Cores without clusters, but with IRAS embedded 
sources had Av NT /Av T = l-Ol;^- Cores with clusters had higher non-thermal linewidths 
still. 

We note that the most complete observations to date suggest that even in low-mass 
cores, the level of turbulent support is comparable to that of thermal support, although 
there are large variations both within and across individual star- forming regions. 1 Moreover, 
there is little evidence that the non-thermal linewidths diminish during the process of star 
formation; if anything, the non-thermal linewidths in cores with IRAS embedded sources 
is higher than in their non-IRAS counterparts. However, we note that the non-thermal 

1 Although these cores are often termed "quiescent", such a designation is somewhat 
misleading. 
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linewidths measured in cores with embedded objects will in general also include infall 
motion towards one or more protostars in a complex flow field. Such infall motions are 
extremely difficult, if not impossible, to separate from turbulent motions on the scale of 
the core without a full spatial mapping of the velocity field, and so it is remains unclear 
(from an observational standpoint) what the timescale for the dissipation of turbulence is, 
or indeed, even whether turbulence decays or increases during star formation. 

2.2. Binary Star Properties 

2.2.1. Binarity Fraction 

Prompted by inconsistencies in multiplicity fractions in earlier surveys (Kuiper 1935, 
1942; Jaschek & Jaschek 1957; Petrie 1960; Jaschek & Gomez 1970), Abt & Levy (1976) 
determined the multiplicity fraction of a sample of 135 bright field stars of types F3-G2, 
and found a binarity fraction (the fraction of all primaries having a detectable secondary 
companion) of 57%. Later, Duquennoy & Mayor (1991) (DM91 hereafter) using a complete 
sample of 164 primary G dwarf stars with types F7-G9, found a remarkably similar binarity 
fraction of 58% (although their conclusions regarding the period distribution differed 
substantially; see below). 

2.2.2. Period Distribution 

DM91 demonstrated that the binary periods in a complete sample of field G dwarf stars 
fits a broad log-normal distribution remarkably well, with a median period of logP^ = 4.8 
and a standard deviation of log0> d = 2.3, where Pd is the period measured in days. DM91's 
result is substantially different than the median period log-P^ = 3.7 found previously by 
Abt & Levy (1976), whose survey was biased by their use of a magnitude-limited sample. 
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Similar period distributions were found for field K dwarfs (Mayor, Duquennoy, Halbwachs, 
& Mermilliod 1992) and field M dwarfs (Fischer & Marcy 1992), though the K dwarf sample 
was found to have a somewhat narrower distribution (log a Pi = 1.9). 

Mathieu (1994) compiled data of pre-main sequence (PMS) binaries extant at that 
time. While the statistics of PMS binaries are not as complete as those in the field, the 
compiled PMS binary period distribution appears to fit a log-normal distribution similar to 
that of Duquennoy & Mayor (1991). More recently, other observers have begun to produce 
tentative evidence in favor of the hypothesis that the period distribution itself varies from 
region to region. For instance, in observations of Upper Scorpius A and B, Brandner 
& Koehler (1998) argue that the period distribution median varies significantly between 
the two regions, with wider binaries being preferred in the low-mass star forming region. 
Similarly, Scally, Clarke, & McCaughrean (1999) argue that a deficiency of wide (a > 1000 
AU) binaries exists in Orion. However, such studies of individual regions and distribution 
tails reduce the counting statistics by necessity (some 20 systems in the case of Brandner 
and Koehler's study; some 3 systems in the case of Scally et al's study), and so we must 
consider them with some caution. 

To date, the period distribution has not been quantitatively explained by theory. 
Mouschovias (1977) envisioned dense cores as non-turbulent, highly subcritical, slowly 
contracting regions initially endowed with angular velocity induced by galactic shear. In 
his picture, magnetic braking becomes inefficient at some critical density at which the core 
becomes supercritical, and thereafter the angular momentum of the gas is conserved during 
the collapse. He invoked variations in this critical density parameter across all star-forming 
regions to explain the broad variation in the binary period distribution. However, it 
remains unclear whether this picture is quantitatively consistent with the observed period 
distribution. 
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More recently, Larson (2002) offered the possibility that tidal interactions among 
binaries in a stellar cluster would tend to transfer angular momentum among the 
binaries, and thereby broaden the period distribution. However, Kroupa & Burkert 

(2001) demonstrated that even under extreme conditions (stellar densities ~ 10 6 pc~ 3 ), 
gravitational encounters were not able to explain the observed breadth of the binary 
period distribution. Indeed, even when a very small amount of broadening of the period 
distribution was obtained in their densest configurations, disruption dominated over 
broadening, so that the binarity fraction dropped from an initial value of unity to ~ .2 
in about 3 • 10 4 yr, in sharp conflict with observation. Indeed, Kroupa & Burkert (2001) 
concluded that the binary period distribution was a remnant of the initial conditions of star 
formation, and was not due to subsequent evolution. However, Bate, Bonnell, & Bromm 

(2002) found that close systems could form in a gas-rich environment through dynamical 
friction. A corollary of Bate's findings is that in a turbulent environment, the angular 
momentum dissipation rate should be sensitive to initial conditions - as had been suggested 
earlier by Larson (2001). We will return to this point later (§6). 

2.2.3. Eccentricity Distribution 

Another important diagnostic of binary stars is the eccentricity distribution. DM91 
divided their binary catalog into short (P < 10 3 d) and long (P > 10 3 d) samples. They 
found that while the short-period binaries had a distribution which peaked at e = .3, 
long-period binaries followed a distribution which was roughly linear (/(e) = 2e), up to 
eccentricities of about .8, at which point the eccentricity distribution tapers off significantly. 

Mathieu (1994) reviewed the data extant at that time for MS and PMS binaries, 
and compared plots of the eccentricity e versus period P for each population of binaries 
separately. One of his principal findings was that the eccentricities of longer period systems 
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(both MS and PMS) were positively correlated with the period of the system, to periods 
of at least 10 3 days. Moreover, the MS eccentricity distribution appears to be roughly 
consistent with the PMS eccentricity distribution. This suggests that the eccentricities of 
MS binaries are established by age of the PMS stage (about 10 6 years) or earlier. 

2.3. Mass Ratio 

An accurate observational determination of the the mass ratio distribution requires a 
number of corrections due to observational biases. DM91 found a mass ratio distribution 
which was consistent with drawing both individual members of the binary system 
independently from the IMF, and could be fitted to a Gaussian. In particular, their mass 
ratio distribution exhibited no maximum towards q — 1. However, a later paper (Mazeh, 
Goldberg, Duquennoy, & Mayor 1992), pointed out that once additional selection-effect 
corrections for spectroscopic binaries with periods less than 3 ■ 10 3 d are taken into account, 
the resultant period distribution for short period binaries is nearly flat, with a slight rise 
towards q — 1. 

3. Initial Molecular Cloud Core Model Formulation 

3.1. Overview 

Our model for turbulent molecular cloud cores will consist of two essential elements. 
The first, the background density model (§3.2) is based upon a Bonnor-Ebert density 
profile of isothermal gas, with pressure support enhanced through turbulent pressure. As 
we will show, this density profile is consistent both with Larson's mean density relationship 
and with resolved observations of molecular cloud cores (§2.1.1). The second element, a 
turbulent velocity field (§3.3), is superposed upon this background density. This velocity 



11 



field is consistent with Larson's linewidth-size relation (§2.1.2). 

In contrast to our assumptions of isothermal gasdynamics and a turbulent velocity 
field, other authors have have assumed that turbulence can effectively be accounted for by 
using a logotropic (McLaughlin & Pudritz 1997) or polytropic pressure component (McKee 
& Holliman 1999), and derived density profiles for their models under the assumption that 
gravity was balanced by the thermal and turbulent pressure of the core. Our approach 
differs fundamentally in that it is essentially a non- equilibrium model - while the kinetic, 
thermal, and gravitational energies are in virial balance, our cores are never static. This 
difference is absolutely critical in the context of the current work, since we seek to tie the 
binary properties to the angular momentum of the turbulent model cores. 



Since our models include turbulent support, we define an effective soundspeed 
c e fj, which includes contributions from both thermal and turbulent pressure, added in 
quadrature. Because pressure depends on only the ID rms of the velocity, we include the 
one-dimensional turbulent support through the effective soundspeed : 



It is both convenient and physically meaningful to scale the radius r to a dimensionless 
unit £ = r/r , where r proportional to the Jeans length at the central density p c and to 
scale the density itself to a dimensionless unit 6 in terms of the central density : 



3.2. Background Density Model 




(2) 




(3) 



ro = 



£ = r/r , 



(4) 
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and 



= p/p, 



•ci 



(5) 



where p e d ge is the edge density of the core, and x — Pel Pedge is the density contrast between 
the center of the core and its edge. Combining the Poisson equation with the momentum 
equation, assuming the pressure support is derived from a polytropic equation of state, and 
demanding both spherical symmetry and force balance in the radial direction, gives us the 
Lane-Emden equation : 



where 7 = 1 + 1/n, and 7 is the ratio of specific heats. We take the boundary conditions 
at £ = as 0(£ = 0) = 1 and ©'(£ = 0) = 0, appropriate for centrally condensed initial 
conditions. While there is no closed- form solution to (6) for 7 = 1, we utilize a simple 
analytic approximation originally suggested by Natarajan & Lynden-Bell (1997): 



where the parameters A, B, C, and D are numerical constants parameterizing the solution. 
Natarajan & Lynden-Bell (1997) suggest A = 50, B = 48, C = VlO, D = ^\2. We assume 
the density distribution is critically stable, with a contrast x — 14 (Bonnor 1957; Ebert 
1955). The resultant approximation accurately describes the critical Bonnor- Ebert sphere 
to within 1% at all spatial points. 

In the absence of turbulence, a critically stable core will have the familiar Bonnor- Ebert 
mass and radius (Bonnor 1957; Ebert 1955) : 




(6) 




(7) 



M th = 1.18 



\/G w p^, 



(8) 



(9) 



-13- 



Where p e d ge is the edge density of the core. 

Using the effective soundspeed in place of the isothermal soundspeed, we have 
expressions for the turbulently-supported Bonnor-Ebert mass and radius: 

c 3 

M BE = 1.18 eff : (10) 

V C 3 Pedge 



R BE = .485 ^1— (11) 

V fledge 



Larson's mean-density law states that the mean density within a given volume scales 
as the inverse size of that region (Larson 1981). Hence the column density through a 
given region should remain roughly independent of the size of the region. Applying this 
mean-density scaling to our Mbe and -Rbe, we find 



M B E /Pcdgc , /10 x 

^ °esd -77- = const. (12) 

BE V 

Hence, because the effective soundspeed c c s must scale with the turbulent linewidth A4, we 
find that the edge density of the cloud core cannot remain a constant, but must scale as 

p * = vi rmm (13) 

where p is the central density, or equivalently, the edge density as one approaches M. — > 0. 
In the remainder of the paper, we will set p = 5 ■ 10~ 20 gm cm" 3 , and c iso = 2 • 10 4 cm 
s _1 , which imply i?BE = -07 pc and Mbe = 2M for a transonic (M. = 1) molecular cloud 
core, which is typical for cores in the Taurus region (Jijina, Myers, & Adams 1999). We 
note that in a large-scale, isothermal, supersonic molecular cloud, shocks will compress and 
rarefy the mean flow of the gas, producing a log-normal distribution in the density field, 
where the width of the distribution depends on the Mach number of the region (Padoan & 
Nordlund 2002). In the current work, we treat only the median density, and not the tails 



14 



of the density distribution. As a result, the minimum mass binary system treated will be 
proportional to the thermal Jeans mass (see §4). This approximation is justified for systems 
whose masses are each close to the median stellar mass, but above the substellar limit. 
We will extend the semi-empirical method to systems containing at least one substellar 
component in a future paper (Fisher 2003). 

3.3. Velocity Field Model 

As we noted earlier (§2.1.2), Larson's linewidth-size relation, as determined by recent 
authors (Jijina, Myers, & Adams 1999) states that the velocity dispersion Av along a line 
of sight through a region of size R scales roughly as R}l 2 on molecular cloud core scales and 



We assume that each Fourier mode in the velocity decomposition is uncorrelated. 
By the Central Limit Theorem, uncorrelated Fourier modes are equivalent to demanding 
that the probability distribution of the amplitude of each mode is drawn from a 
Gaussian distribution, and the phase of each mode is drawn from a uniform distribution. 
Decomposing an individual mode 5k in Fourier space in terms of its amplitude and phase as 
Sk = rkexp(i(pk), the Gaussian probability distribution gk for the amplitude r\ and phase 
<pk can be written as 



The standard deviation Ok is referred to as the power spectrum of the perturbation 
spectrum. As is conventional, we write P(k) — cr k . By examining the successive moments 
of a Gaussian distribution, one can show that only the first two moments - the mean and 
standard deviation - are needed to entirely specify the distribution. Hence, specifying 
the mean and standard deviation completely determines the Gaussian spectrum. In the 



above. 




(14) 
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following discussion, we assume that the mean of the perturbing field is zero, and that the 
perturbing field is homogeneous and isotropic, so that the amplitude of the power spectrum 
depends solely on the magnitude of the A;— vector. 

By Parseval's Theorem, the total power introduced in /c-space is identical to the total 
power in real space. Hence, the 3D rms velocity Av in a 3D spherical volume in real space 
extending from L min to L max in the radial direction can be easily evaluated by integrating 
the power spectrum in k-space from k m i n = 27r/L max to A; max = 27r/L min : 



Adopting a power-law power spectrum of the form P(k) = Ak ™, then integrating over 
angles and wavenumbers, we have 



Assuming that most of the power is on large scales (n > 3) and that L min <C L max , the 
maximum wavenumber k max corresponding to the minimum scale in the problem can be 
neglected, so we have 



When n = 4, we can recover Larson's Law (Av oc Lmax)- This n = 4 Gaussian velocity 
perturbation power spectrum provides the basis for our turbulent models. Specifically, we 
assume all background velocities are zero, and apply three separate turbulent velocity fields 
to our model molecular cloud core via 




(15) 




(16) 




(17) 



v x (r) 



5 x (r)c iso M 



(18) 



v y {r) 



v z (r) 



5 y (r)c iso M 
$z(r)c iso M 



(19) 



(20) 
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where S x (r), S y (r), and 5 z (r) are three realizations of a A; 4 Gaussian field, Cj SO is the 
isothermal soundspeed, and M. = Av/c iso is the 3D rms turbulent Mach number. 

While Gaussian perturbation fields and the structures which arise from them have 
been a topic of active investigation in the cosmological context for several decades, our use 
of them in the context of star formation differs in two critical respects. First, application 
of a perturbation to the density field requires that the total power must be small : a 
nonlinear perturbation with sufficient total power will drive the density negative in some 
regions. By applying perturbations to the velocity field directly, the total power can be 
made arbitrarily large without encountering unphysical negative densities. Second, whereas 
cosmological density fluctuations typically introduce most power on small scales (the 
inflationary Harrison-Zel'dovich spectrum uses n — 1), we introduce most power on large 
scales, in accord with observations of star forming regions. 

One should appropriately ask how the ratio of rotational to gravitational binding 
energy, (3, depends on the 3D turbulent Mach number Ai, while keeping the realization of 
the turbulence fixed. Interestingly enough, under the assumption of a critical Bonnor-Ebert 
density distribution, for transonic and supersonic cores, a fixed realization yields a (5 that is 
completely independent of the edge-density of the core, and only weakly dependent on the 
turbulent Mach number. Specifically, 

P«( J 4)(^) (2D 



I J \M 2 
M 2 

(5 oc —s (23) 

ie, for transonic and supersonic cores with Ai > 1, (3 varies by roughly a factor of 3 from 
Ai = 1 to Ai = oo. Here I is the moment of inertia of the initial core. This scaling 
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explains why turbulent core models naturally produce the same median value of (3 as is 
seen in observation without any fine-tuning of parameters: in critical Bonnor-Ebert spheres 
dominated by turbulent pressure support, Mach scaling applies, and as a result the models 
are scale-free. Therefore, in the supersonic regime, for a fixed turbulent realization, (3 
does not depend on any model parameters. This scaling holds approximately even down 
to the transonic regime, so that the value of (3 is primarily determined from the slope 
of the turbulent spectrum. The robust agreement between the predicted value of (3 and 
observation gives strong support for our model of the turbulent velocity field. 

Similar Gaussian turbulent spectra imposed on the velocity field have been used in 
a variety of numerical simulations of turbulence in the interstellar medium. Dubinski, 
Narayan, & Phillips (1995) were apparently the first to suggest that Larsonian turbulence 
in the interstellar medium could be generated using Gaussian random velocity fields with 
index n = 4. A number of authors, beginning with Gammie & Ostriker (1996), studied 
the evolution of turbulence in molecular clouds using the same turbulent spectrum - 
initially establishing an Alfvenic spectrum of waves by an initial velocity fluctuation on a 
constant density and magnetic field background. Burkert & Bodenheimer (2000) computed 
linewidth gradients induced by turbulent eddies on the scale of molecular cloud cores. 
Hujeirat, Myers, Camenzind, & Burkert (2000) computed the initial decay of turbulence 
and subsequent collapse (followed to an early evolutionary time) of molecular cloud cores 
supported by Alfvenic turbulence. Lastly, Bate, Bonnell, & Bromm (2002) computed the 
evolution of a large, highly supersonic (M ~ 7) turbulent core. 

4. Semi-Empirical Model for Binary Formation via Turbulent Fragmention 

We may now predict the distribution of binary periods implied by our isolated turbulent 
core models. Since our knowledge of the detailed physics of turbulent fragmentation is still 
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rudimentary at best, in order to predict binary properties from our initial molecular cloud 
core models, we must make several semi-empirical assumptions motivated by observation. 
We will assume that the star-formation efficiency e* = (Ml + M 2 )/M core is constant. 
Further, in direct analogy to the star-formation efficiency, we define the conversion efficiency 
of angular momentum of cores to binaries, ej = J / J cor e, and assume it is constant as well. 
In addition, we will assume that the individual stellar masses are uncorrelated, so that both 
may be randomly drawn from the IMF. Lastly, we will assume that the orbital eccentricities 
of all binaries are drawn from a distribution /(e) = 2e. 

We recognize that these assumptions may not be rigorously correct; in particular, 
based on observations of field stars, there is evidence (Duquennoy & Mayor 1991) that short 
(P < 10 3 d) period binary systems follow a different eccentricity distribution than long 
period binaries, and may also have a different mass-ratio distribution (Mazeh, Goldberg, 
Duquennoy, & Mayor 1992). However, because such short-period systems will be strongly 
affected by disk-star interactions, it is likely that their orbital parameters undergo significant 
evolution. Hence it, in the context of obtaining initial binary properties in the current work, 
we will explore the consequences of the hypothesis that all binary systems are formed with 
uncorrelated masses drawn from the IMF, and obeying a thermal eccentricity distribution. 

We proceed as follows. 

1) First, we randomly draw two masses Mi and M 2 from the IMF. We uniformly draw 
£ over the interval from 0-1 twice, and assign masses (scaled to solar) and mass ratio q 
according to 

ry £72 I ry £74 

-(0 = -08 + 7l ^_ + J,f 8 , (24) 

where 71 = .19, 72 = 1.55, 73 = .050, 74 = .6 (Kroupa, Gilmore, & Tout 1991). We note 
that this IMF introduces a cutoff at .O8M , and so does not extend into the substellar 
range. We define Mi to be the greater of these two masses; hence q = M 2 /Mi. 
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2) Next, using the star formation efficiency factor, set the mass of the initial protostellar 

core: 



Since we now know the parent core mass as well as the thermal Bonnor-Ebert mass, we 
may also set the 3D turbulent Mach number, using equation (3.2) : 



In cases where M core < M t h, the computed core mass is less than the thermal Bonnor-Ebert 
value, implying that the initial core is inconsistent with the selected binary masses and 
assumed star formation efficiency In these instances, we simply reject the model binary 
and draw both stars again. 

3) We then stochastically generate three perturbation cubes, and perturb the velocity 
field of our model core using the known Mach number Ai, thereby setting the resultant 
core angular momentum J corc as well as the binary system angular momentum J = ejJ core . 

4) Next, we draw the eccentricity of the binary orbit from the thermal distribution 
/(e) = 2e. We do this by drawing a number £ from a uniform distribution ranging from 
to 1, and setting e = y 7 ^. 

5) Finally, knowing J, e, M and q, we may compute the period P and semi-major axis 
a of the binary system, from 



(Mi + M 2 ) 



(25) 



core 




(26) 




(27) 



and 




(28) 



(See Appendix A for derivation.) 
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These relations are often written down for the special case of an equal-mass, circular 
binary (e = 0, and q — 1). Clearly, in that case, the larger the value of J/M, the greater the 
amount of angular momentum in the system, and hence, the longer the period. In addition, 
the smaller the total mass M of the system, the less the influence of gravity, and hence, 
the wider the binary. However, although it is commonly not recognized, for fixed mass and 
angular momentum, the eccentricity and the mass ratio can also have a substantial impact 
on the orbital properties of binaries. For instance, for a fixed angular momentum, the 
larger the eccentricity e of a system, the wider the binary. Similarly, for a fixed mass, the 
smaller the mass ratio q, the wider the binary. Although apparently trivial, these scaling 
relationships will have an important bearing on the statistical properties of binary star 
systems, as we shall soon see. 

Our models introduce two free parameters: e* and ej, the mass and angular momentum 
star formation efficiencies, respectively. The first parameter is relatively well-constrained by 
both theory and observation (Matzner & McKee 2000). For a given e*, the ej parameter 
is obtained by requiring that the resultant model binary period distribution median agree 
with the observed value. In effect, the imposition of this constraint reduces the parameter 
space to one free parameter, which we take to be e*. We note, from equation (27), that with 

5 /3 

the constraint that the model period median must agree with observation, that ej oc e* . 
The nonlinearity of the dependence of ej on e* has a simple interpretation: for low star 
formation efficiency, a larger amount of mass must be accreted onto the binary, which 
implies a greater angular momentum loss efficiency due to the larger turbulent linewidth on 
larger scales. 

Note that our model implicitly assumes that the mass of the system M, the mass 
ratio q, and the orbital eccentricity e are uncorrelated, whereas the angular momentum 
J (and hence the period P) is explicitly dependent upon M : the higher the mass of 
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the binary system, the greater the mass and Mach number M. of the initial core. These 
assumptions differ from those adopted by Kroupa (1995) in his formulation of the inverse 
binary population synthesis problem, where he assumed that P, q, and e were uncorrelated. 
In contrast, in our formulation, the period P depends explicitly on q and e through angular 
momentum conservation (equation 27), which implies that our model predicts non-trivial 
correlations between P and q, and between P and e. We return to these points in later 
sections (§5.2 and §5.3), where we detail the nature of these correlations. 

5. Results 

5.1. Binary Period Distribution 

Figure 1 shows the binary period distribution obtained for the star formation efficiency 
e* = 0.26. The standard deviation of the computed distribution is <jp d = 1.7; slightly 
narrower than that of the field, though consistent with the PMS distribution. The computed 
distribution is in remarkably good agreement with Mathieu's cataloged visual binaries with 
periods > 10 4 d. The computed distribution exhibits a somewhat smaller binarity fraction 
at periods shorter than 10 3 d than either the field or the PMS distribution, although given 
the numbers involved, this result is of marginal statistical significance. 

There is a small, but systematic increase in the width of the distribution with increasing 
star formation efficiency. We interpret this increasing width as a consequence of the decrease 
in angular momentum loss efficiency with increasing star formation efficiency. As a result, 
higher star formation efficiency models have slightly wider distributions of pre-stellar core 
specific angular momenta. The distribution of initial core specific angular angular momenta 
for the case of e* = 0.26 is shown in figure 2. 
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5.2. Eccentricity Versus Period 

As we noted in §4, in our semi-empirical model of turbulent star formation, the angular 
momentum of the binary system is directly related to the initial angular momentum of 
the turbulent core. As a result, the period and eccentricity are no longer uncorrelated 
quantities, and the period is explicitly dependent upon the eccentricity; for a fixed angular 
momentum and mass, the greater the eccentricity, the longer the period. From observation 
(§2.2.3), we know that in binaries wider than the tidal circularization limit (P ~ 10 d), the 
period tends to be correlated with the eccentricity, which is in fact a robust feature of our 
model (see equation 27 ). 

We plot eccentricity versus the log (base 10) of the period for in figure 5, for the case 
of intermediate star-formation efficiency (e* = 0.5). For comparison, binaries from the field 
(Duquennoy & Mayor 1991), and from PMS regions (Mathieu 1994), are shown in star and 
plus symbols, respectively, in the plot on the right. The vertical dashed line indicates the 
tidal circularization limit found by Duquennoy & Mayor (1991) at about P d ~ 11 d. The 
key result here is that the model systems exhibit a positive correlation between eccentricity 
and period, which is qualitatively similar to that observed. However, it is certainly true 
that our models predict an overabundance of highly eccentric binaries at shorter periods. 
It is quite likely that additional angular momentum dissipation mechanisms, including 
coupling to circumstellar disks, may help circularize such systems. 

5.3. Mass Ratio Versus Period 

The semi-empirical model also predicts a correlation between mass ratio and period. 
The reason for the correlation is, once again, simply related to conservation of angular 
momentum. A system with two unequal masses, and a small mass ratio (q <^ 1) will have 
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a longer period than a similar system with the same total angular momentum and mass, 
but with equal masses. The results of 200 model systems for the case of intermediate 
star-formation efficiency (e* = 0.5) are shown in figure 4. 

We emphasize that in order to explain the correlation between mass ratio and 
period, the semi-empirical model framework does not require that the masses themselves be 
correlated. Instead, the model suggests that it is possible to form two stars independently 
within a single turbulent core, with the resultant period anti-correlated with the mass ratio 
simply through angular momentum conservation. 

5.4. Binary System Mass Versus Semimajor Axis 

As we noted previously, observations of low-mass systems have revealed that such 
systems do not follow the same statistical trends observed in solar-mass and higher binaries. 
Perhaps the most noticeable difference is the trend for low-mass systems to occur in 
tighter binaries. In the semi-empirical model, two factors contribute to the explanation 
of this effect. First, lower- mass binary systems originate from parent cores of lower mass, 
and hence, lower specific angular momentum. Second, drawing stars independently from 
the IMF, a low-mass system will naturally tend to have a more equal-mass ratio. Taken 
together, the model predicts that low-mass systems should naturally tend to occur in 
shorter-period systems, which is in fact what is observed (Martin & Basri 2001). 

In figures 6 - 7, we display the results from two model systems of variable star formation 
efficiency. The minimum system mass represented for a given star-formation efficiency 
e* is indicated with a dashed horizontal line; no model systems form below the line. For 
comparison, we also plot the largest observable binary separations a maX) as determined by 
Close, Siegler, Freed, Biller (2003) in their figure 15. They found that the two straight lines 
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of constant orbital velocity a max = 23.2(M tot /.185) AU and a max = 1000(M tot /.185) AU 
(where M tot is in solar masses) describe the maximal separation upper envelopes relatively 
well. 

We note that, above the horizontal minimum mass line, our model systems fit the 
Close et al. upper envelopes remarkably well; very few systems lie to the right of the 
envelopes, quite independent of the star-formation efficiency In the context of the current 
paper, our assumption of a fixed Jeans mass requires that low mass systems must form in 
regions of low star-formation efficiency Hence, it appears that a combination of a variety 
of star-formation efficiencies can explain the results of Close et al. More general models 
including constant star-formation efficiency and supersonic turbulence on large scales may 
also be able to explain the observations; a topic we are investigating further (Fisher 2003). 



6. Discussion 

One limitation of our current approach is the choice of a fixed central density po- By 
specifying p , we have implicitly restricted our attention to systems above a mass e*M th . 
As we pointed out previously (§3.2), shock compressions in a turbulent molecular cloud 
will produce a log-normal distribution of po- We expect that a set of binary systems 
formed from cores with varying M t h resulting from to shock compression will in fact have a 
slightly broader period distribution than that computed here. Similarly, the IMF adopted 
introduces a cutoff at .O8M , even though substellar objects are also formed in star-forming 
regions. These effects may partially account for the slightly smaller period distribution 
standard deviations computed here (<7p d ~ 1.6 — 2.1), in comparison to that found in the 
field (a Pd ~ 2.3) (Duquennoy & Mayor 1991). 

Detailed theoretical work treating wind-driven outflows (Matzner & McKee 2000) 
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predicts that the star-formation efficiency on the molecular cloud scale should lie in the 
range 25 - 70%. For star-formation efficiencies in this range, our models predict core 
properties consistent with observation. In particular, with our fiducial scaling, the median 
Mach number M lies in the range of 1.2 and 0.9 over the range e* = 0.1 — 0.9 - quite close 
to the actual median value of .9 observed in Taurus (Jijina, Myers, & Adams 1999). We 
note that since both our pre-stellar core and final binary properties are relatively insensitive 
to the star-formation efficiency (the median prestellar core Mach number M. varies only by 
a factor of 2 for star-formation efficiencies over the range e* = 0.1 — 0.9) our models do not 
provide further constraints on the star-formation efficiency itself. 

A key consequence of the semi-empirical model for star formation described here stems 
from the fact that many observational properties can be quantitatively described with a 
constant angular momentum loss factor. In contrast, several authors have advocated a 
"chaotic" description star formation (Larson 2001; Bate, Bonnell, & Bromm 2002), in which 
a variety of cluster interactions play a key role in setting the initial properties of binaries. 
In this paper, we have explicitly demonstrated that chaotic interactions are not required 
to explain many properties of binaries wider than P > 10 3 d, provided that turbulent 
fragmentation can in fact directly produce the mass and eccentricity distributions which we 
have assumed. Note that we distinguish the fact that our initial conditions are turbulent, 
and therefore inherently stochastic (though deterministic) from the sensitivity upon initial 
conditions implied by a chaotic model. We would, however, agree that dynamical friction 
is likely to play an important role in the formation of tight binaries, and that subsequent 
tidal interactions may also act to broaden the initial period distribution found here. 

Another key consequence of our model is the importance which angular momentum loss 
from initial molecular cloud cores plays in determining the properties of binaries. Although 
this problem has been identified by a number of authors (Simon 1992; Bodenheimer 1995; 
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Larson 2002), virtually all simulations of multiple star formation on the molecular cloud 
core scale done to date remain purely hydrodynamic, and have neglected the influence of 
the magnetic field entirely. 2 As a result, the only means of angular momentum transport 
treated in most calculations is gravitational torquing. Given the importance of angular 
momentum loss in determining binary properties, and the relative effectiveness of magnetic 
braking in axisymmetric simulations of single-star formation (Basu & Mouschovias 1994), 
caution must be applied in comparing the results of purely hydrodynamic simulations 
directly against observation; in all likelihood, the magnetic field plays a critical role in 
determining the binary angular momentum. 

RTF would like to thank Richard Klein and Christopher McKee for their prescient 
suggestion to investigate the physical role of turbulence in multiple star formation as a 
thesis topic. Thanks also go to Marc Davis and Matt Craig for the use of their Gaussian 
perturbation generation program. This research has made use of NASA's Astrophysics 
Data System Bibliographic Service. 

A. Characteristic Scales 

In this appendix, we derive the characteristic scales of a binary system consisting of 
two masses Mi and M 2 , in terms of its angular momentum J, its mass ratio q = M 2 /Mi, 
and its orbital eccentricity e. 

From standard two-body classical mechanics in an inverse square law potential 

2 To the best of our knowledge, the only exceptions are the work of Hujeirat, Myers, 
Camenzind, & Burkert (2000) and Boss (2000), although Boss did not include a fully self- 
consistent treatment of the MHD equations. 
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(Goldstein 1980), we know the angular momentum of a binary system is completely 
specified by the semi-major and semi-minor axes of the bound orbit, as well as the reduced 
mass of the two-body system \i = M 1 M 2 /(M 1 + M 2 ), and the inverse square law constant 
of proportionality k : 



a 



J= 7T72V^ ( A1 ) 



where \x = M 1 M 2 /(M 1 + M 2 ) is the reduced mass of the two body system, and k = GM 1 M 2 
for the gravitational force. Hence, since b = ay/1 — e 2 , we can determine the final orbital 
separation as a function of the eccentricity e, total mass M = Mi +M 2 , and specific angular 
momentum J/M : 

Defining the ratio of secondary to primary mass q = M 2 /Mi, for e < 1, the semi-major axis 
of a bound orbit is 

_1/J\ 2 1 1 (1 + g) 4 
a ~G\M) l^M^q^~ [A6) 

For a typical molecular cloud core specific angular momentum value of J/M = 3- 10 20 cm 2 /s, 



3- 10 20 cm 2 / S y 1-e 2 \ M J q 
Then, from Kepler's third law, we can determine the period of the final binary system : 

n 27r 1 ( J V (! + <1? 1 



G 2( 1 _ e 2)3/2 ^ M J g 3 M 2 



J/M \ 3 (1 + g) 6 fM Q x :: 



(A5) 



P = 300 ^rs ir -^^ 31020 cm2/ J -^^J (A6) 

We anticipate that explicit dependence of P on e in equation (A6) may cause some readers, 
who are accustomed to the fact that the Keplerian period depends only on a, and not on e, 
some undue consternation. These concerns will be rapidly alleviated when it is recognized 
that, for fixed masses and angular momentum, a more eccentric binary system must have 
a wider separation (eqn. A4). Our equation (A6) expresses the period in terms of the 
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angular momentum J in favor of the semimajor axis a, and hence the origin of the explicit 
dependence of P on e. 

Lastly, using the classical two-body inverse square law result for the total energy E of 
the system 

E = -± (A7) 
la 

one can also determine the final energy E of the binary system : 

E = k*W (A8 » 
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e* 


ej 


logP d 




e*M th /M e 


a/AU 


j (cm 2 s 1 ) 


M 


0.1 


.002 


4.9 


1.6 


.14 


30 


1.3 • 10 21 


1.7 


0.3 


.01 


4.9 


1.7 


.31 


34 


7.7 • 10 20 


1.2 


0.5 


.04 


4.8 


1.8 


.51 


30 


4.5 • 10 20 


1.0 


0.7 


.06 


4.9 


2.0 


.72 


35 


3.5 • 10 20 


.89 


0.9 


.10 


5.0 


2.1 


.92 


45 


3.4 • 10 20 


.87 







Table 1: A table of model parameters, for a variety of star formation efficiencies (e* = 
0.1 — 0.9) fitted such that the median model period agrees with that of the field. Shown are 
the mass and angular momentum star- formation efficiencies e* and ej, the median of the log 
of the period distribution (in days) \ogP d , the standard deviation of the log of the period 
distribution (in days) logop d , the minimum mass binary system (in solar masses) e^M t h/M Q) 
the median of the semi-major axis (in AU) a, the median prestellar core angular momentum 
j (in cm 2 s _1 ), and the median prestellar core Mach number M. 
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Fig. 1. — Histogram of logP^. The left panel shows our numerical results for 200 model 
systems for the case of star-formation efficiency e* = 0.26. For comparison, the right-hand 
panel shows the period distribution inferred from PMS stars (Mathieu 1994) and from field 
stars (Duquennoy & Mayor 1991), in solid and dashed lines, respectively. 
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Fig. 2. — Histogram of the log of the model initial core specific angular mom entum j = J/M, 
measured in cm 2 s" 1 , for 200 binary-producing cores, shown for a low star- format ion efficiency 
(e, = 0.26). 
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Fig. 3. — Histogram of loga(AU). The numerical results for 200 model systems for the cases 
of low star formation efficiency (e* = 0.26), are shown, with Poisson error bars drawn. 
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Fig. 4. — Binary mass ratio versus log Pd- The circles show the results of 200 model systems 
for the cases of intermediate star-formation efficiency (e* = 0.5). 
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Fig. 5. — Eccentricity versus log Pa- The circles show the results of 200 model systems for the 
case of intermediate star-formation efficiency (e* = 0.5). For comparison, binaries from the 
field (Duquennoy & Mayor 1991), and from PMS regions (Mathieu 1994), are shown in star 
and plus symbols, respectively, in the plot on the right. The vertical dashed line indicates 
the tidal circularization limit found by Duquennoy & Mayor (1991) at about Pa ~ 11 d. 
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Fig. 6. — Log semimajor axis versus log total mass. The circles show the results of model 
systems for the very low star- formation effficiency case e* = 0.1. The two solid slanted 
lines indicate the upper envlopes indicated by Close, Siegler, Freed, Biller (2003) (see text). 
The horizontal dashed line is drawn to delineate the minimum mass binary system for this 
star-formation efficiency 
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Fig. 7. — Log semimajor axis versus log total mass, as in fig. 6. The circles show the results 
of model systems for the intermediate star formation effficiency case e* = 0.5. 



